#!/usr/bin/python
# -*- coding: ascii -*-
#
from __future__ import print_function, division
import os, sys
from os import path
import numpy as np
import CoolProp
import glob
from warnings import warn
from time import clock

import CoolProp.constants
from CoolProp.CoolProp import PropsSI,generate_update_pair,get_parameter_index,set_debug_level,get_phase_index
from CoolProp import AbstractState as State

import matplotlib as mpl
import matplotlib.pyplot as plt
import matplotlib.patches as mpatches
import copy
from itertools import cycle
from matplotlib import gridspec, ticker
#from jopy.dataPlotters import roundList, range_brace

def range_brace(x_min, x_max, mid=0.5,
                beta1=50.0, beta2=100.0, height=1,
                initial_divisions=11, resolution_factor=1.5):
    """
    http://stackoverflow.com/questions/1289681/drawing-braces-with-pyx
    x,y = range_brace(0, 100)
    ax.plot(x, y,'-')
    ax.plot(y, x,'-')
    """
    # determine x0 adaptively values using second derivitive
    # could be replaced with less snazzy:
    #   x0 = NP.arange(0, 0.5, .001)
    x0 = np.array(())
    tmpx = np.linspace(0, 0.5, initial_divisions)
    tmp = beta1**2 * (np.exp(beta1*tmpx)) * (1-np.exp(beta1*tmpx)) / np.power((1+np.exp(beta1*tmpx)),3)
    tmp += beta2**2 * (np.exp(beta2*(tmpx-0.5))) * (1-np.exp(beta2*(tmpx-0.5))) / np.power((1+np.exp(beta2*(tmpx-0.5))),3)
    for i in range(0, len(tmpx)-1):
        t = int(np.ceil(resolution_factor*max(np.abs(tmp[i:i+2]))/float(initial_divisions)))
        x0 = np.append(x0, np.linspace(tmpx[i],tmpx[i+1],t))
    x0 = np.sort(np.unique(x0)) # sort and remove dups
    # half brace using sum of two logistic functions
    y0 = mid*2*((1/(1.+np.exp(-1*beta1*x0)))-0.5)
    y0 += (1-mid)*2*(1/(1.+np.exp(-1*beta2*(x0-0.5))))
    # concat and scale x
    x = np.concatenate((x0, 1-x0[::-1])) * float((x_max-x_min)) + x_min
    y = np.concatenate((y0, y0[::-1])) * float(height)
    return (x,y)

#try:
    #from jopy.dataPlotters import BasePlotter
    #bp = BasePlotter()
#except:
    #bp = None

bp = None


# The basic settings for he plots
xypoints = 1000
loops    = 1
repeat   = 1
runs     = 0
maxruns  = 5
plot     = True
calc     = True
check    = True
folder   = "dataTTSE"
figures  = "figuresTTSE"
#np.random.seed(1984)


fluids  = ["CO2","Pentane","R134a","Water","Air","LiBr-0%"]
fluids  = ["CO2","Pentane","R134a","Water"]
fluids  = ["Air"]

#glskeys = [r"\glsentryshort{co2}",r"\glsentryshort{pentane}",r"\glsentryshort{r134a}",r"\glsentryshort{water}",r"\glsentryshort{air}",r"\glsentryshort{libr} \SI{0}{\percent}"]
#glskeys = [r"\ce{CO2}",r"n-Ppentane",r"R134a",r"Water",r"Air",r"\glsentryshort{libr} \SI{0}{\percent}"]
repList = []
#for i in range(len(fluids)):
#    repList.append(fluids[i])
#    repList.append(glskeys[i])

#backends = ["INCOMP","HEOS","REFPROP"]
backends = ["HEOS","REFPROP"]
backends = ["HEOS"]

#repList.append("HEOS")
#repList.append(r"\glsentryshort{cp}")
#repList.append("REFPROP")
#repList.append(r"\glsentryshort{rp}")


#CoolProp.CoolProp.set_debug_level(51)

pStr = path.dirname(path.abspath(__file__))
fStr = path.splitext(path.basename(__file__))[0]
def getFolderName():
    folderName = path.join(pStr,folder)
    if not path.isdir(folderName):
        print("Creating data directory "+folderName)
        os.makedirs(folderName)
    return folderName

def getFigureFolder():
    folderName = path.join(pStr,figures)
    if not path.isdir(folderName):
        print("Creating data directory "+folderName)
        os.makedirs(folderName)
    return folderName

repList.append("TimeComp-")
repList.append("chapters/FluidProperties/"+path.basename(getFigureFolder())+"/TimeComp-")

def getFileName(qualifiers=[]):
    fileName = path.join(getFolderName(),"-".join(qualifiers))
    return fileName

# Some file handling
def loadNpzData(backend,fluid):
    dicts = {}
    globber = getFileName([backend,fluid])+'_[0-9][0-9][0-9].npz'
    for fname in glob.glob(globber):
        dataDict = dict(np.load(fname))
        dicts[str(dataDict["name"])] = dataDict
    #if len(dicts)<1:
    #    #print("No readable file found for {0}".format(globber))
    #    dataDict = dict(name=str(0).zfill(3))
    #    dicts[str(dataDict["name"])] = dataDict
    return dicts

def saveNpzData(backend,fluid,dicts,start=0,stop=-1):
    keys = dicts.keys()
    keys.sort()
    for k in keys[start:stop]:
        data = dicts[k]
        fname = getFileName([backend,fluid])+'_{0}.npz'.format(str(data['name']).zfill(3))
        np.savez(fname, **data)
    return True

def splitFluid(propsfluid):
    fld = propsfluid.split("::")
    if len(fld)==2:
        backend =  fld[0]
        fld = fld[1]
    else:
        backend = None
        fld = fld[0]
    fld = fld.split("-")
    if len(fld)==2:
        conc  = float(fld[1].strip('%'))/100.0
        fld = fld[0]
    else:
        conc  = None
        fld = fld[0]
    return backend,fld,conc

def getInpList(backend):
    if backend=="HEOS": return ["DT","HP"]
    elif backend=="REFPROP": return ["DT","HP"]
    elif backend=="INCOMP": return ["PT","HP"]
    else: raise ValueError("Unknown backend.")

def getOutList(inp=None):
    if inp == "HP":
        return [["Tmax"],["D"],["S"],["T"],["D","S","T"]]
    elif inp == "DT":
        return [["Tmax"],["H"],["P"],["S"],["H","P","S"]]
    elif inp == "PT":
        return [["Tmax"],["H"],["D"],["S"],["H","D","S"]]
    else:
        raise ValueError("Unknown inputs.")

#def getPhaseString(iPhase):
#    for i in range(11):
#        if getPhaseString(i)==sPhase:
#            return i
#     if iPhase==1: return "liquid"
#     elif iPhase==2: return "supercritical"
#     elif iPhase==3: return "supercritical_gas"
#     elif iPhase==4: return "supercritical_liquid"
#     elif iPhase==5: return "critical_point"
#     elif iPhase==6: return "gas"
#     elif iPhase==7: return "twophase"
#     elif iPhase==8: return "unknown"
#     elif iPhase==9: return "not_imposed"
#     else: raise ValueError("Couldn't find phase.")

def getPhaseNum(sPhase):
    return get_phase_index(sPhase)
#     for i in range(11):
#         if getPhaseString(i)==sPhase:
#             return i

def getOutKey(out): return "".join(out)

def getOutLabel(out): return ",".join(out)

def getTimeKey(inp,out): return "_".join([inp,getOutKey(out)])

def getVectorKey(inp,out): return getTimeKey(inp, out)+"_V"

def getCriticalProps(propsfluid):
    backend,_,_ = splitFluid(propsfluid)
    if backend!="INCOMP":
        p_crit_m = PropsSI('pcrit',"T",0,"D",0,propsfluid)  *0.995
        T_crit_m = PropsSI('Tcrit',"T",0,"D",0,propsfluid)  *1.005
        d_crit_m = PropsSI('rhocrit',"T",0,"D",0,propsfluid)*0.995
        h_crit_m = PropsSI('H',"T",T_crit_m,"D",d_crit_m,propsfluid)
        s_crit_m = PropsSI('H',"T",T_crit_m,"D",d_crit_m,propsfluid)
    else:
        p_crit_m = None
        T_crit_m = None
        d_crit_m = None
        h_crit_m = None
        s_crit_m = None
    return dict(P=p_crit_m,T=T_crit_m,D=d_crit_m,H=h_crit_m,S=s_crit_m)


def getPTRanges(propsfluid):
    backend,_,_ = splitFluid(propsfluid)
    # Setting the limits for enthalpy and pressure
    T_min = PropsSI('Tmin',"T",0,"D",0,propsfluid)+1
    T_max = PropsSI('Tmax',"T",0,"D",0,propsfluid)-1

    if backend == "REFPROP":
        T_min = max(T_min,PropsSI('Ttriple',"T",0,"D",0,propsfluid))+1
        p_min = PropsSI(  'P',"T",T_min,"Q",0,propsfluid)+1
        p_max = PropsSI(   'pmax',"T",0,"D",0,propsfluid)-1
    elif backend == "INCOMP":
        p_min =   1.5*1e5
        p_max = 200.0*1e5
    else:
        T_min = max(T_min,PropsSI('Ttriple',"T",0,"D",0,propsfluid))+1
        p_min = PropsSI('ptriple',"T",0,"D",0,propsfluid)
        p_min = max(p_min,PropsSI('pmin',"T",0,"D",0,propsfluid))+1
        p_max = PropsSI(   'pmax',"T",0,"D",0,propsfluid)-1

    # One more check to debug things:
    #p_min = max(p_min,0.01e5)
    #T_min = max(T_min,200)
    #p_max = min(p_max,200e5)
    #T_max = min(T_max,1750)

    p_range  = np.logspace(np.log10(p_min),np.log10(p_max),xypoints)
    T_range  = np.linspace(T_min,T_max,xypoints)
    return p_range,T_range

    #p_max = min(PropsSI('pcrit',"T",0,"D",0,fluid)*20, p_max)
    #T_max = min(PropsSI('Tcrit',"T",0,"D",0,fluid)* 3, T_max)


def getLists(propsfluid):
    backend,_,_ = splitFluid(propsfluid)
    """Returns randomised lists of all properties within the ranges"""
    p,T = getPTRanges(propsfluid)
    p_min = np.min(p)
    p_max = np.max(p)
    T_min = np.min(T)
    T_max = np.max(T)

    if backend=="INCOMP":
        h_min = PropsSI('H','T',T_min,'P',p_min,propsfluid)
        h_max = PropsSI('H','T',T_max,'P',p_max,propsfluid)
    else:
        critProps = getCriticalProps(propsfluid)
        h_min = PropsSI('H','T',T_min,'Q',0,propsfluid)
        h_max = PropsSI('H','T',T_min,'Q',1,propsfluid)
        h_max = max(PropsSI('H','T',critProps["T"],'D',critProps["D"],propsfluid),h_max)
        h_max = (h_max-h_min)*2.0+h_min

    loop = True
    count = 0
    while loop:
        count += 1

        h_list = np.random.uniform(          h_min,          h_max,int(xypoints*2.0))
        p_list = np.random.uniform(np.log10(p_min),np.log10(p_max),int(xypoints*2.0))
        p_list = np.power(10,p_list)

        out = ["T","D","S"]
        res = PropsSI(out,"P",p_list,"H",h_list,propsfluid)
        T_list = res[:,0]
        d_list = res[:,1]
        s_list = res[:,2]

        mask = np.isfinite(T_list) & np.isfinite(d_list) & np.isfinite(s_list)
        if np.sum(mask)<xypoints:
            if False:
                print(h_list); print(p_list); print(T_list); print(d_list); print(s_list)
            print("There were not enough valid entries in your result vector: {0:d} > {1:d} - rerunning".format(xypoints,np.sum(mask)))
            loop = True
        else:
            loop = False
            p_list = p_list[mask][0:xypoints]
            h_list = h_list[mask][0:xypoints]
            T_list = T_list[mask][0:xypoints]
            d_list = d_list[mask][0:xypoints]
            s_list = s_list[mask][0:xypoints]
            return dict(P=p_list,T=T_list,D=d_list,H=h_list,S=s_list)

        maxTries = 4
        if count > maxTries:
            loop = False
            raise ValueError("{0}: Could not fill the lists in {0} runs, aborting.".format(propsfluid,maxTries))

def getInpValues(inp,dataDict):
    in1 = inp[0]
    in2 = dataDict[in1]
    in3 = inp[1]
    in4 = dataDict[in3]
    return in1,in2,in3,in4

def getStateObj(propsfluid):
    backend,fld,conc = splitFluid(propsfluid)
    # fluidstr holds the full information and fluid is only the name
    # Initialise the state object
    if backend is not None:
        state = State(backend,fld)
    else:
        state = State(fld)
    #if backend=="INCOMP":
    #    state.set_mass_fractions([0.0])
    if conc is not None:
        try:
            state.set_mass_fractions([conc])
        except:
            pass
    return state

def getSpeedMeas(out,in1,in2,in3,in4,propsfluid,vector=False):
    pair, out1, _ = generate_update_pair(get_parameter_index(in1),in2[0],get_parameter_index(in3),in4[0])
    if out1==in2[0]: swap=False
    else: swap=True
    if swap:
        input1 = in4
        input2 = in2
    else:
        input1 = in2
        input2 = in4
    state = getStateObj(propsfluid)
    outList = [get_parameter_index(i) for i in out]
    resLst = np.empty((repeat,))
    timLst = np.empty((repeat,))
    if vector:
        for j in range(repeat):
            timLst.fill(np.inf)
            lrange = range(len(input1))
            resTmp = np.inf
            if len(outList)==1 and outList[0]==CoolProp.constants.iT_max:
                t1=clock()
                for l in lrange:
                    for o in outList:
                        resTmp = state.keyed_output(o)
                t2=clock()
                timLst[j] = (t2-t1)*1e6/float(len(input1))
            else: # We have to update before doing other things
                t1=clock()
                for l in lrange:
                    state.update(pair,input1[l],input2[l])
                    for o in outList:
                        resTmp = state.keyed_output(o)
                t2=clock()
                timLst[j] = (t2-t1)*1e6/float(len(input1))
        res = None
        tim = np.min(timLst) # Best of (repeat)
        return res,tim
    else:
        res = np.empty_like(input1)
        res.fill(np.inf)
        tim = np.empty_like(input1)
        tim.fill(np.inf)
        for i in range(len(input1)):
            resLst.fill(np.inf)
            timLst.fill(np.inf)
            for j in range(repeat):
                lrange = range(loops)
                resTmp = np.inf
                if len(outList)==1 and outList[0]==CoolProp.constants.iT_max:
                    t1=clock()
                    for _ in lrange:
                        for o in outList:
                            resTmp = state.keyed_output(o)
                    t2=clock()
                    timLst[j] = (t2-t1)*1e6/float(loops)
                    resLst[j] = resTmp
                else: # We have to update before doing other things
                    inV1 = input1[i]
                    inV2 = input2[i]#*(1.0+(l/1000.0)*pow(-1,l)) for l in lrange ]
                    t1=clock()
                    for l in lrange:
                        state.update(pair,inV1,inV2)
                        for o in outList:
                            resTmp = state.keyed_output(o)
                    t2=clock()
                    timLst[j] = (t2-t1)*1e6/float(loops)
                    resLst[j] = resTmp
            if not np.all(resLst==resLst[0]):
                raise ValueError("Not all results were the same.")
            res[i] = resLst[0]
            tim[i] = np.min(timLst) # Best of three (repeat)
        return res,tim


def checkDataSet(propsfluid,dataDict,fill=True,quiet=False):
    if not check: return
    backend,_,_ = splitFluid(propsfluid)
    if not quiet: print("\n\n-- {0:16s} --".format(propsfluid),end="")
    # Test for required inputs
    newInputs = False
    inLists = getLists(propsfluid)
    for inp in getInpList(backend):
        if not quiet: print("\n{0:2s}: ".format(inp),end="")
        for inVal in inp:
            if inVal not in dataDict: # A problem
                if not fill:
                    raise ValueError("The input {0:1s} is missing or faulty, cannot continue.".format(inVal))
                #dataDict[inVal] = inLists[inVal]
                dataDict.update(inLists)
                newInputs = True
                if not quiet: print("{0:s}*({1:d}),".format(inVal,len(dataDict[inVal])),end="")
            else:
                if not quiet: print("{0:s} ({1:d}),".format(inVal,len(dataDict[inVal])),end="")
        # All inputs are there
        in1,in2,in3,in4 = getInpValues(inp,dataDict)
        #in2 = in2[:3]
        #in4 = in4[:3]
        if in2.shape != in4.shape:
            raise ValueError("The stored data for {0:s} and {1:s} do not have the same shape.".format(in1,in3))
        if in2.shape != inLists[in1].shape:
            raise ValueError("The stored data for {0:s} and its list do not have the same shape {1} vs {2}.".format(in1,in2.shape,inLists[in1].shape))
        # Check for time data
        for out in getOutList(inp):
            key = getTimeKey(inp, out)
            okey = getOutKey(out)
            if key not in dataDict or newInputs or not np.all(np.isfinite(dataDict[key])):
                if not fill:
                    raise ValueError("The time data for {0:s} is missing or faulty, cannot continue.".format(key))
                res,tim = getSpeedMeas(out, in1, in2, in3, in4, propsfluid)
                dataDict[key] = tim
                dataDict[okey] = res # We calculated in, why not use it here...
                if not quiet: print("{0:s}*({1:d}),".format(key,len(dataDict[key])),end="")
            else:
                if not quiet: print("{0:s} ({1:d}),".format(key,len(dataDict[key])),end="")
            if dataDict[key].shape != in2.shape or not np.all(np.isfinite(dataDict[key])):
                raise ValueError("The stored time data for {0:s} does not have the same shape as the inputs.".format(key))
        # Check for vectors
        for out in getOutList(inp):
            key = getVectorKey(inp, out)
            if key not in dataDict or not np.all(np.isfinite(dataDict[key])):
                if not fill:
                    raise ValueError("The fluid data for {0:s} is missing or faulty, cannot continue.".format(key))
                res,tim = getSpeedMeas(out, in1, in2, in3, in4, propsfluid, vector=True)
                dataDict[key] = tim
                if not quiet: print("{0:s}*({1:d}),".format(key,dataDict[key].size),end="")
            else:
                if not quiet: print("{0:s} ({1:d}),".format(key,dataDict[key].size),end="")
            if dataDict[key].size!=1 or not np.all(np.isfinite(dataDict[key])):
                raise ValueError("The vector data for {0:s} does not have the correct size {1}..".format(key,dataDict[key].size))

#     inp = getInpList(backend)[0] # Explicit calls
#      # Check for properties
#      for out in getOutList(inp)[:-1]:
#          key = getOutKey(out)
#          if key not in dataDict or not np.all(np.isfinite(dataDict[key])):
#              if not fill:
#                  raise ValueError("The fluid data for {0:s} is missing or faulty, cannot continue.".format(key))
#              res = PropsSI(out,in1,in2,in3,in4,propsfluid)
#              dataDict[key] = res
#              if not quiet: print("{0:s}*({1:d}),".format(key,len(dataDict[key])),end="")
#          else:
#              if not quiet: print("{0:s} ({1:d}),".format(key,len(dataDict[key])),end="")
#          if dataDict[key].shape != in2.shape or not np.all(np.isfinite(dataDict[key])):
#              raise ValueError("The stored data for {0:s} does not have the same shape as the inputs {1} vs {2}..".format(key,dataDict[key].shape,in2.shape))
#     # Check for phase
#     for out in ["Phase"]:
#         if backend!="HEOS":
#             dataDict[key] = np.zeros_like(a, dtype, order, subok)
#         key = getOutKey(out)
#         if key not in dataDict or newInputs or not np.all(np.isfinite(dataDict[key])):
#             if not fill:
#                 raise ValueError("The phase data for {0:s} is missing or faulty, cannot continue.".format(key))
#             res = np.empty_like(in2)
#             res.fill(np.inf)
#             for i in range(len(in2)):
#                 res[i] = PropsSI(out,in1,in2[i],in3,in4[i],propsfluid)
#             dataDict[key] = res
#             if not quiet: print("{0:s}*({1:d}),".format(key,len(dataDict[key])),end="")
#         else:
#             if not quiet: print("{0:s} ({1:d}),".format(key,len(dataDict[key])),end="")
#         if dataDict[key].shape != in2.shape or not np.all(np.isfinite(dataDict[key])):
#             raise ValueError("The stored data for {0:s} does not have the same shape as the inputs {1} vs {2}..".format(key,dataDict[key].shape,in2.shape))
#
#             # Now we use the vector data
#             key = getVectorKey(inp, out)
#             if key not in dataDict or not np.all(np.isfinite(dataDict[key])):
#                 if not fill:
#                     raise ValueError("The vector data for {0:s} is missing or faulty, cannot continue.".format(key))
#                 dataDict[key] = np.empty_like(in2)
#                 dataDict[key].fill(np.inf)
#                 res = []
#                 for _ in range(repeat):
#                     t1=clock()
#                     PropsSI(out,in1,in2,in3,in4,propsfluid)
#                     t2=clock()
#                     res.append((t2-t1)/float(len(in2)))
#                 dataDict[key] = np.min(res)
#                 if not quiet: print("{0:s}*({1}),".format(key,dataDict[key]),end="")
#             else:
#                 if not quiet: print("{0:s} ({1}),".format(key,dataDict[key]),end="")
#             try:
#                 float(dataDict[key])
#             except:
#                 raise ValueError("The stored vector data for {0:s} cannot be casted to float.".format(key))
#        if not quiet: print("")

# All data is loaded and checked, we can calculate more now
def getEntryCount(dicts,backend,fld):
    return len(fluidData[fld][backend].keys())

def getUKey(fld,bck,inp,out):
    return "-".join([fld,bck,inp,"".join(out)])

def getData(fld,backend,inp,out,fluidData):
    inputs1 = []
    inputs2 = []
    values  = []
    times   = []
    i1key = inp[0]
    i2key = inp[1]
    vkey  = getOutKey(out)
    tkey  = getTimeKey(inp,out)
    for dkey in fluidData[fld][backend]:
        cData = fluidData[fld][backend][dkey]
        inputs1.append(cData[i1key])
        inputs2.append(cData[i2key])
        values.append( cData[vkey] )
        times.append(  cData[tkey] )
    ret = {}
    if len(inputs1)>0:
        ret[i1key] = np.concatenate(inputs1)
        ret[i2key] = np.concatenate(inputs2)
        ret[vkey]  = np.concatenate(values )
        ret[tkey]  = np.concatenate(times  )
    return ret

def getSingleData(fld,backend,key,fluidData):
    #print("Getting: "+fld+", "+backend+", "+key)
    values  = []
    for dkey in fluidData[fld][backend]:
        if key in fluidData[fld][backend][dkey]:
            if "P" in fluidData[fld][backend][dkey]:
                # TODO: Fix this, do we need the mask?
                #mask = fluidData[fld][backend][dkey]["P"]>0.3e5
                mask = fluidData[fld][backend][dkey]["P"]>0.0e5
                try:
                    values.append( fluidData[fld][backend][dkey][key][mask] )
                except Exception as e:
                    values.append( fluidData[fld][backend][dkey][key] )
                    print(e)
                    pass
            else:
                values.append( fluidData[fld][backend][dkey][key] )
    if len(values)>0:
        if np.size(values[0])>1:
            return np.concatenate(values)
        else:
            return np.array(values)
    return None

def fillDict(fld,backend,fluidData,curDict,curKeys):
    if curDict is None: curDict = {}
    for key in curKeys:
        vals = getSingleData(fld, backend, key, fluidData)
        if vals is not None: curDict[key] = vals
    return curDict

################################################

fluidData = {}

for fluidstr in fluids:
    _,fld,_ = splitFluid(fluidstr)
    if fld not in fluidData: fluidData[fld] = {}
    for backend in backends:
        if backend not in fluidData[fld]: # Try to add it
            propsfluid = "::".join([backend,fluidstr])
            try:
                PropsSI('Tmax',"T",0,"D",0,propsfluid)
                fluidData[fld][backend] = loadNpzData(backend,fld)
            except:
                pass
            if backend in fluidData[fld]:
                for k in fluidData[fld][backend]:
                    checkDataSet(propsfluid, fluidData[fld][backend][k], fill=False, quiet=True)
        else: # Already added backend for fluid
            pass

lastSave = 0
while runs<maxruns and calc:
    check = True # force checking for new records
    runs +=1
    # Now we have the data structure with the precalculated data
    for fluidstr in fluids:
        _,fld,_ = splitFluid(fluidstr)
        for backend in fluidData[fld]:
            propsfluid = "::".join([backend,fluidstr])
            dicts      = fluidData[fld][backend]
            keys = list(dicts.keys())
            keys.sort()
            while len(keys)<runs:
                if len(keys)<1: newKey = 0
                else: newKey = int(keys[-1])+1
                if newKey>999:
                    raise ValueError("Too many dicts: {0}>999".format(newKey))
                k = str(newKey).zfill(3)
                dicts[k] = {}
                dicts[k]['name']=k
                try:
                    checkDataSet(propsfluid, dicts[k], fill=True, quiet=False)
                except Exception as e:
                    print("There was an error, dataset {0} from {1}.might be faulty:\n{2}".format(k,propsfluid,str(e)))
                    pass
                todel = []
                #for k in dicts:
                try:
                    checkDataSet(propsfluid, dicts[k], fill=False, quiet=True)
                except Exception as e:
                    print("There was an error, removing dataset {0} from {1}.\n{2}".format(k,propsfluid,str(e)))
                    todel.append(k)
                for t in todel: del dicts[t]
                keys = list(dicts.keys())
                keys.sort()

    # Updated all dicts for this backend, saving data
    if runs>=maxruns or (lastSave+4)<runs:
        print("\n\nDone processing fluids, saving data: ")
        for fld in fluidData:
            for backend in fluidData[fld]:
                saveNpzData(backend, fld, fluidData[fld][backend], start=lastSave, stop=runs)
                print("{0} ({1})".format(backend+"::"+fld,len(fluidData[fld][backend].keys()[lastSave:runs])),end=", ")
        print("")
        lastSave = runs


if not plot: sys.exit(0)

# exclLst = [["Tmax"],["H"],["D"],["S"],["H","D","S"],["D"],["S"],["T"],["D","S","T"]]
#
### Start with a temporary dictionary that holds all the data we need
# for fld in fluidData:
#     cstData = {} # Data from calls to constants (overhead)
#     expData = {} # Data from explicit EOS calls
#     impData = {} # Data from EOS calls that require iterations
#     for backend in fluidData[fld]:
#
#         curKeys = []
#         for inp in getInpList(backend):
#             for out in getOutList(inp)[:1]:
#                 curKeys.append(getTimeKey(  inp, out))
#                 curKeys.append(getVectorKey(inp, out))
#         curDict = {}
#         fillDict(fld,backend,fluidData,curDict,curKeys)
#         cstData[backend] = curDict
#
#         curKeys = []
#         for inp in getInpList(backend)[:1]:
#             for out in getOutList(inp)[1:]:
#                 curKeys.append(getTimeKey(  inp, out))
#                 curKeys.append(getVectorKey(inp, out))
#         curDict = {}
#         fillDict(fld,backend,fluidData,curDict,curKeys)
#         expData[backend] = curDict
#
#         curKeys = []
#         for inp in getInpList(backend)[1:]:
#             for out in getOutList(inp)[1:]:
#                 curKeys.append(getTimeKey(  inp, out))
#                 curKeys.append(getVectorKey(inp, out))
#         curDict = {}
#         fillDict(fld,backend,fluidData,curDict,curKeys)
#         impData[backend] = curDict


#         curDict[backend] = {}
#         for key in :
#             vals = getSingleData(fld, backend, key, fluidData)
#             if vals is not None: curDict[backend][key] = vals
#             if curDict
#
#         cstData[backend] = {}
#         for out in getOutList(inp[0])[0]:
#             res = getData(fld,backend,inp[0],out,fluidData)
#             cstData[backend].update(res)
#         for out in getOutList(inp[1])[0]:
#             res = getData(fld,backend,inp[1],out,fluidData)
#             cstData[backend].update
#
#         expData[backend] = {}
#         for out in getOutList(inp[0])[1:]:
#             res = getData(fld,backend,inp[0],out,fluidData)
#             expData[backend].update(res)
#
#         impData[backend] = {}
#         for out in getOutList(inp[1])[1:]:
#             res = getData(fld,backend,inp[1],out,fluidData)
#             impData[backend].update(res)

    #############################################################
    # All data is available in the dicts now.
    #############################################################
    # The first thuing to do is to print some statistical
    # measures to give you an idea about the data.
#     try:
#         #dataOHCP = [cstData["HEOS"]["DT_Tmax"]   , cstData["HEOS"]["HP_Tmax"]   ]
#         #dataOHRP = [cstData["REFPROP"]["DT_Tmax"], cstData["REFPROP"]["HP_Tmax"]]
#         print("\n{0} - {1} points ".format(fld,np.size(cstData["HEOS"]["DT_Tmax"])))
#         print("Overhead CoolProp: {0:5.3f} us".format(np.mean(cstData["HEOS"]["DT_Tmax"])))#,   np.mean(cstData["HEOS"]["HP_Tmax"]))
#         print("Overhead REFPROP : {0:5.3f} us".format(np.mean(cstData["REFPROP"]["DT_Tmax"])))#,np.mean(cstData["REFPROP"]["HP_Tmax"]))
#         print("Mean EOS in  CoolProp: f(rho,T): {0:9.3f} us, f(h,p): {1:9.3f} us".format(np.mean(expData["HEOS"]["DT_HPS"]),np.mean(impData["HEOS"]["HP_DST"])))
#         print("Std. dev. in CoolProp: f(rho,T): {0:9.3f} us, f(h,p): {1:9.3f} us".format(np.std(expData["HEOS"]["DT_HPS"]) ,np.std(impData["HEOS"]["HP_DST"])))
#         print("Minimum in   CoolProp: f(rho,T): {0:9.3f} us, f(h,p): {1:9.3f} us".format(np.min(expData["HEOS"]["DT_HPS"]) ,np.min(impData["HEOS"]["HP_DST"])))
#         print("Maximum in   CoolProp: f(rho,T): {0:9.3f} us, f(h,p): {1:9.3f} us".format(np.max(expData["HEOS"]["DT_HPS"]) ,np.max(impData["HEOS"]["HP_DST"])))
#         print("Mean EOS in   REFPROP: f(rho,T): {0:9.3f} us, f(h,p): {1:9.3f} us".format(np.mean(expData["REFPROP"]["DT_HPS"]),np.mean(impData["REFPROP"]["HP_DST"])))
#         print("Std. dev. in  REFPROP: f(rho,T): {0:9.3f} us, f(h,p): {1:9.3f} us".format(np.std(expData["REFPROP"]["DT_HPS"]) ,np.std(impData["REFPROP"]["HP_DST"])))
#         print("Minimum in    REFPROP: f(rho,T): {0:9.3f} us, f(h,p): {1:9.3f} us".format(np.min(expData["REFPROP"]["DT_HPS"]) ,np.min(impData["REFPROP"]["HP_DST"])))
#         print("Maximum in    REFPROP: f(rho,T): {0:9.3f} us, f(h,p): {1:9.3f} us".format(np.max(expData["REFPROP"]["DT_HPS"]) ,np.max(impData["REFPROP"]["HP_DST"])))
#         print("")
#
#         print("\n{0} - {1} points ".format(fld,np.size(cstData["HEOS"]["DT_Tmax_V"])))
#         print("Overhead CoolProp: {0:5.3f} us".format(np.mean(cstData["HEOS"]["DT_Tmax_V"])))#,   np.mean(cstData["HEOS"]["HP_Tmax"]))
#         print("Overhead REFPROP : {0:5.3f} us".format(np.mean(cstData["REFPROP"]["DT_Tmax_V"])))#,np.mean(cstData["REFPROP"]["HP_Tmax"]))
#         print("Mean EOS in  CoolProp: f(rho,T): {0:9.3f} us, f(h,p): {1:9.3f} us".format(np.mean(expData["HEOS"]["DT_HPS_V"]),np.mean(impData["HEOS"]["HP_DST_V"])))
#         print("Std. dev. in CoolProp: f(rho,T): {0:9.3f} us, f(h,p): {1:9.3f} us".format(np.std(expData["HEOS"]["DT_HPS_V"]) ,np.std(impData["HEOS"]["HP_DST_V"])))
#         print("Minimum in   CoolProp: f(rho,T): {0:9.3f} us, f(h,p): {1:9.3f} us".format(np.min(expData["HEOS"]["DT_HPS_V"]) ,np.min(impData["HEOS"]["HP_DST_V"])))
#         print("Maximum in   CoolProp: f(rho,T): {0:9.3f} us, f(h,p): {1:9.3f} us".format(np.max(expData["HEOS"]["DT_HPS_V"]) ,np.max(impData["HEOS"]["HP_DST_V"])))
#         print("Mean EOS in   REFPROP: f(rho,T): {0:9.3f} us, f(h,p): {1:9.3f} us".format(np.mean(expData["REFPROP"]["DT_HPS_V"]),np.mean(impData["REFPROP"]["HP_DST_V"])))
#         print("Std. dev. in  REFPROP: f(rho,T): {0:9.3f} us, f(h,p): {1:9.3f} us".format(np.std(expData["REFPROP"]["DT_HPS_V"]) ,np.std(impData["REFPROP"]["HP_DST_V"])))
#         print("Minimum in    REFPROP: f(rho,T): {0:9.3f} us, f(h,p): {1:9.3f} us".format(np.min(expData["REFPROP"]["DT_HPS_V"]) ,np.min(impData["REFPROP"]["HP_DST_V"])))
#         print("Maximum in    REFPROP: f(rho,T): {0:9.3f} us, f(h,p): {1:9.3f} us".format(np.max(expData["REFPROP"]["DT_HPS_V"]) ,np.max(impData["REFPROP"]["HP_DST_V"])))
#         print("")
#
#     except Exception as e:
#         print(str(e.message))
#         pass
#
# try:
#     fld = 'Water'
#     cstData = {} # Data from calls to constants (overhead)
#     expData = {} # Data from explicit EOS calls
#     impData = {} # Data from EOS calls that require iterations
#     for backend in fluidData[fld]:
#         curKeys = []
#         for inp in getInpList(backend):
#             for out in getOutList(inp)[:1]:
#                 curKeys.append(getTimeKey(  inp, out))
#                 curKeys.append(getVectorKey(inp, out))
#         curDict = {}
#         fillDict(fld,backend,fluidData,curDict,curKeys)
#         cstData[backend] = curDict
#
#         curKeys = []
#         for inp in getInpList(backend)[:1]:
#             for out in getOutList(inp)[1:]:
#                 curKeys.append(getTimeKey(  inp, out))
#                 curKeys.append(getVectorKey(inp, out))
#         curDict = {}
#         fillDict(fld,backend,fluidData,curDict,curKeys)
#         expData[backend] = curDict
#
#         curKeys = []
#         for inp in getInpList(backend)[1:]:
#             for out in getOutList(inp)[1:]:
#                 curKeys.append(getTimeKey(  inp, out))
#                 curKeys.append(getVectorKey(inp, out))
#         curDict = {}
#         fillDict(fld,backend,fluidData,curDict,curKeys)
#         impData[backend] = curDict
#
#
#     print("Done")

def autolabel(ax,rects,means,stds,lens=0,fac=1):
    return
    # attach some text labels
    yerr = (stds*100.0)/means
    ypos = np.max(means) #+ np.max(stds)
    for i in range(len(rects)):
        xpos = rects[i].get_x()+rects[i].get_width()/2.
        #ax.text(xpos, 1.05*ypos, '{0:s}{1:4.2f}{2:s}{3:3.1f}{4:s}'.format(r'',means[i]*fac,r'us +- ',yerr[i],r'%'), rotation=90, ha='center', va='bottom', fontsize='smaller')
        ax.text(xpos, 1.05*ypos, '{0:s}{1:4.2f}{2:s}{3:3.1f}{4:s}'.format(r'\SI{',means[i]*fac,r'}{\us} (\SI{\pm',yerr[i],r'}{\percent})'), rotation=90, ha='center', va='bottom', fontsize='smaller')
        #ax.text(xpos, 0.25*ypos, str(lens[i]), rotation=90, ha='center', va='bottom', fontsize='smaller')


def axislabel(txt,ax,xPos,yPos=-1):
    ax.text(xPos, yPos, txt, rotation=45, ha='center', va='top', fontsize='xx-small')


#############################################################
# All data is available in the dicts now.
#############################################################
# The first plot contains the time data, this is averaged and
# plotted as a bar graph with the standard deviation.
hatchLst    = ["", "///"  ,"\\\\\\" ]
backendsLst = ["INCOMP","HEOS","REFPROP"]
for fluidstr in fluids[:-1]:
    _,fld,_ = splitFluid(fluidstr)
    outCst = []
    labCst = []
    outExp = []
    labExp = []
    outImp = []
    labImp = []
    DEBUG=True
    for backend in backendsLst:
        # Backend exists in fluid data?
        try:
            for i,inp in enumerate(getInpList(backend)):
                for j,out in enumerate(getOutList(inp)):
                    if j==0: # First output is Tmax
                        if backend not in fluidData[fld]:
                            outCst.append([0])
                            labCst.append("Dummy")
                            if DEBUG: print("Added a dummy for {0} and {1},{2},{3}".format(fld,backend,inp,out))
                        else:
                            outCst.append(getSingleData(fld, backend, getTimeKey(inp, out), fluidData))
                            labCst.append(getTimeKey(inp, out))
                        continue
                    elif i==0: # First input is explicit
                        if backend not in fluidData[fld]:
                            outExp.append([0])
                            labExp.append("Dummy")
                            if DEBUG: print("Added a dummy for {0} and {1},{2},{3}".format(fld,backend,inp,out))
                        else:
                            outExp.append(getSingleData(fld, backend, getTimeKey(inp, out), fluidData))
                            labExp.append(getTimeKey(inp, out))
                        continue
                    elif i==1:
                        if backend not in fluidData[fld]:
                            outImp.append([0])
                            labImp.append("Dummy")
                            if DEBUG: print("Added a dummy for {0} and {1},{2},{3}".format(fld,backend,inp,out))
                        else:
                            outImp.append(getSingleData(fld, backend, getTimeKey(inp, out), fluidData))
                            labImp.append(getTimeKey(inp, out))
                        continue
                    else:
                        raise ValueError("Wrong number.")
        except Exception as e:
            print(e)
            sys.exit(1)

    # Do the plotting
    if bp is not None:
        bp.figure = None
        fg = bp.getFigure()
        ccycle = bp.getColorCycle(length=3)
    else:
        fg = plt.figure()
        ccycle = cycle(["b","g","r"])

    #fg = plt.figure()
    ax1 = fg.add_subplot(111)
    ax2 = ax1.twinx()
    rects1  = []
    labels1 = []
    rects2  = []
    labels2 = []
    rects3  = []
    labels3 = []
    col1 = ccycle.next()
    col2 = ccycle.next()
    col3 = ccycle.next()


    numBackends = len(backendsLst)

    step   = 1
    width  = step/(numBackends+1)  # the width of the bars
    offset = -0.5*numBackends*width

    entries = int((len(outCst)+len(outExp)+len(outImp))/numBackends)

#     for o in range(entries):
#         ids = np.empty((numBackends,))
#         for b in range(numBackends):
#             i = o*numBackends+b
#             ids[b] = offset + o*step + b*width
#             j = i - 0
#             if j < len(outCst):
#                 rects1.extend(ax1.bar(ids[b], np.mean(outCst[j]), width, color=col1, hatch=hatchLst[b]))#, yerr=np.std(curList[i]), ecolor='k'))
#
#
#             j = i - 0
#             if j < len(outCst):
#                 rects1.extend(ax1.bar(ids[b], np.mean(outCst[j]), width, color=col1, hatch=hatchLst[b]))#, yerr=np.std(curList[i]), ecolor='k'))
#             else:
#                 j = i-len(outCst)
#                 if j < len(outExp):
#                     rects2.extend(ax1.bar(ids[b], np.mean(outExp[j]), width, color=col2, hatch=hatchLst[b]))#, yerr=np.std(curList[i]), ecolor='k'))
#                 else:
#                     j = i-len(outCst)-len(outExp)
#                     if j < len(outImp):
#                         rects3.extend(ax2.bar(ids[b], np.mean(outImp[j]), width, color=col3, hatch=hatchLst[b]))#, yerr=np.std(curList[i]), ecolor='k'))
#                     else:
#                         raise ValueError("Do not go here!")
    DEBUG=True

    entries = 2
    for o in range(entries):
        ids = np.empty((numBackends,))
        for b in range(numBackends):
            i = b*entries+o
            try:
                ids[b] = offset + o*step + b*width
                rects1.extend(ax1.bar(ids[b], np.mean(outCst[i]), width, color=col1, hatch=hatchLst[b], rasterized=False))#, yerr=np.std(curList[i]), ecolor='k'))
                if DEBUG:
                    print("Plotting {0}: {1:7.1f} - {2} - {3}".format(fld,np.mean(outCst[i]),"cst.",b))
                    #print(ids[b],labCst[i])
                    #axislabel(labCst[i], ax1, ids[b]+0.5*width)
            except:
                pass

    offset += entries*step
    #entries = int(len(outExp)/numBackends)
    entries = 4
    for o in range(entries):
        ids = np.empty((numBackends,))
        for b in range(numBackends):
            i = b*entries+o
            try:
                ids[b] = offset + o*step + b*width
                rects2.extend(ax1.bar(ids[b], np.mean(outExp[i]), width, color=col2, hatch=hatchLst[b], rasterized=False))#, yerr=np.std(curList[i]), ecolor='k'))
                if DEBUG:
                    print("Plotting {0}: {1:7.1f} - {2} - {3}".format(fld,np.mean(outExp[i]),"exp.",b))
                    #print(ids[b],labExp[i])
                    #axislabel(labExp[i], ax1, ids[b]+0.5*width)
            except:
                pass

    x_newaxis = np.max(ids)+1.5*width
    plt.axvline(x_newaxis, color='k', linestyle='dashed')

    offset += entries*step
    entries = 4
    for o in range(entries):
        ids = np.empty((numBackends,))
        for b in range(numBackends):
            i = b*entries+o
            try:
                ids[b] = offset + o*step + b*width
                rects3.extend(ax2.bar(ids[b], np.mean(outImp[i]), width, color=col3, hatch=hatchLst[b], rasterized=False))#, yerr=np.std(curList[i]), ecolor='k'))
                if DEBUG:
                    print("Plotting {0}: {1:7.1f} - {2} - {3}".format(fld,np.mean(outImp[i]),"imp.",b))
                    #print(ids[b],labImp[i])
                    #axislabel(labImp[i], ax1, ids[b]+0.5*width)
            except:
                pass


    #ax1.set_xlim([ids.min()-2.5*width,ids.max()+2.5*width])
    ax1.spines['top'].set_visible(False)
    ax2.spines['top'].set_visible(False)
    ax1.xaxis.set_ticks_position('bottom')
    ax2.xaxis.set_ticks_position('bottom')

    labels=[r"ex.",r"im.",r"$h$",r"$\rho|p$",r"$s$",r"all",r"$\rho$",r"$s$",r"$T$",r"all"]
    ax1.set_xticks(range(len(labels)))
    ax1.set_xticklabels(labels)
    #ax1.yaxis.get_label().set_verticalalignment("baseline")

    x_min = rects1[0].get_x()
    dx = rects1[0].get_width()
    x_max = rects3[-1].get_x()

    x_min = x_min - 1*dx
    x_max = x_max + 2*dx

    ax1.set_xlim([x_min,x_max])

    y_min = 0

    y_max_c = np.nanmax([a.get_height() for a in rects1])
    y_max_e = np.nanmax([a.get_height() for a in rects2])
    y_max_i = np.nanmax([a.get_height() for a in rects3])

    y_max = np.max([y_max_c,y_max_e,y_max_i/10.0])

    y_max = np.ceil(1.3*y_max/10.0)*10.0

    ax1.set_ylim([y_min,y_max])
    ax2.set_ylim([y_min,y_max*10.0])

    ratio = 10.0/4.0 * y_max / 250.0 # height of 10 for 4 points if y_max==250

    x_min = rects1[ 0].get_x()
    x_max = rects1[-1].get_x() + dx
    x,y = range_brace(x_min, x_max)
    dy  = np.ceil(y_max_c/10.0)*10.0
    y   = dy+y*ratio*(x[-1]-x[0])
    ax1.plot(x, y, ls='-',color='k')
    ax1.text(np.mean(x), np.max(y), "const.", rotation=0, ha='center', va='bottom', fontsize='medium')

    x_min = rects2[ 0].get_x()
    x_max = rects2[-1].get_x() + dx
    x,y = range_brace(x_min, x_max)
    dy  = np.ceil(y_max_e/10.0)*10.0
    y   = dy+y*ratio*(x[-1]-x[0])
    ax1.plot(x, y, ls='-',color='k')
    ax1.text(np.mean(x), np.max(y), "explicit", rotation=0, ha='center', va='bottom', fontsize='medium')

    x_min = rects3[ 0].get_x()
    x_max = rects3[-1].get_x() + dx
    x,y = range_brace(x_min, x_max)
    dy  = np.ceil(y_max_i/100.0)*10
    y   = dy+y*ratio*(x[-1]-x[0])
    ax1.plot(x, y, ls='-',color='k')
    ax1.text(np.mean(x), np.max(y), "implicit", rotation=0, ha='center', va='bottom', fontsize='medium')

    #ax1.text(x_newaxis*0.9, y_max*0.9, "<- left axis", rotation=0, ha='right', va='bottom', fontsize='medium')
    #ax1.text(x_newaxis*1.1, y_max*0.9, "right axis ->", rotation=0, ha='left', va='bottom', fontsize='medium')

    handles = []
    for h in (rects1[0], rects2[1], rects3[2]):
        handles.append(copy.copy(h))
        handles[-1].set_facecolor('white')
        handles.append(copy.copy(h))
        handles[-1].set_hatch('')
    labels = (r'$p,T$-fit', r'constant', r'CoolProp', r'explicit, $f(p|\rho,T)$', r'REFPROP', r'implicit, $f(h,p)$')

    if bp is not None:
        bp.drawLegend(ax=ax1,
          loc='lower center',
          bbox_to_anchor=(0.5, 1.05),
          ncol=3,
          handles=handles,
          labels=labels)
    else:
        ax1.legend(handles,labels,
          loc='lower center',
          bbox_to_anchor=(0.5, 1.05),
          ncol=3)


    ax1.set_ylabel(r'Time per explicit call (us)')
    ax2.set_ylabel(r'Time per implicit call (us)')
    fg.savefig(path.join(getFigureFolder(),"TimeComp-"+fld.lower()+".pdf"))
    if bp is not None:
        ax1.set_ylabel(r'Time per explicit call (\si{\us})')
        ax2.set_ylabel(r'Time per implicit call (\si{\us})')
        mpl.rcParams['text.usetex'] = True
        fg.savefig(path.join(getFigureFolder(),"TimeComp-"+fld.lower()+"-tex.pdf"))
        mpl.rcParams['text.usetex'] = False
        # Fix the wrong baseline
        for tick in ax1.get_xaxis().get_major_ticks():
            tick.set_pad(2*tick.get_pad())
            tick.label1 = tick._get_text1()
        for lab in ax1.xaxis.get_ticklabels():
            lab.set_verticalalignment("baseline")
            #lab.set_pad(1.5*lab.get_pad())


        #ax1.set_xticklabels(labels)
        #
        #for tick in ax1.xaxis.get_major_ticks():
        #    tick.label1.set_horizontalalignment('center')
        bp.savepgf(path.join(getFigureFolder(),"TimeComp-"+fld.lower()+".pgf"),fg,repList)

    plt.close()

# for fld in fluids:
#     try:
#         if bp is not None:
#             bp.figure = None
#             fg = bp.getFigure()
#             ccycle = bp.getColorCycle(length=3)
#         else:
#             fg = plt.figure()
#             ccycle = cycle(["b","g","r"])
#
#         #fg = plt.figure()
#         ax1 = fg.add_subplot(111)
#         ax2 = ax1.twinx()
#
#         if "INCOMP" in fluidData[fld]:
#             el = 3
#             hatch = ["","//","x"]
#         else:
#             el = 2
#             hatch = ["//","x"]
#
#         #one standard deviation above and below the mean of the data
#         width  = 0.25        # the width of the bars
#         step   = 1
#         offset = -step-0.5*el*width
#
#         lab   = []
#         rects1 = []
#         rects2 = []
#
#         curList = []
#         if el==3: curList.append(getSingleData(fld, "INCOMP" , "PT_Tmax", fluidData)*1e2)
#         curList.append(getSingleData(fld, "HEOS"   , "DT_Tmax", fluidData)*1e2)
#         curList.append(getSingleData(fld, "REFPROP", "DT_Tmax", fluidData)*1e2)
#
#         lab.extend(["Tmax"])
#
#         offset += step
#         curN   = len(curList)
#         curInd = [ offset + i * width for i in range(curN)]
#         curCol = ccycle.next()
#         for i in range(len(hatch)):
#             rects1.extend(ax1.bar(curInd[i], np.mean(curList[i]), width, color=curCol, hatch=hatch[i]))#, yerr=np.std(curList[i]), ecolor='k'))
#         autolabel(ax1,rects1[0:],np.mean(curList,axis=1),np.std(curList,axis=1),fac=1e-2)
#
#         curList = []
#         if el==3: curList.append(getSingleData(fld, "INCOMP" , "PT_H", fluidData))
#         curList.append(getSingleData(fld, "HEOS"   , "DT_H", fluidData))
#         curList.append(getSingleData(fld, "REFPROP", "DT_H", fluidData))
#         lab.extend(["H"])
#         offset += step
#         curN   = len(curList)
#         curInd = [ offset + i * width for i in range(curN)]
#         curCol = ccycle.next()
#         for i in range(len(hatch)):
#             rects1.extend(ax1.bar(curInd[i], np.mean(curList[i]), width, color=curCol, hatch=hatch[i]))#, yerr=np.std(curList[i]), ecolor='k'))
#         autolabel(ax1,rects1[el:],np.mean(curList,axis=1),np.std(curList,axis=1))
#
#         curList = []
#         if el==3: curList.append(getSingleData(fld, "INCOMP" , "PT_D", fluidData))
#         curList.append(getSingleData(fld, "HEOS"   , "DT_P", fluidData))
#         curList.append(getSingleData(fld, "REFPROP", "DT_P", fluidData))
#         if el==3: lab.extend(["D/P"])
#         else: lab.extend(["P"])
#         offset += step
#         curN   = len(curList)
#         curInd = [ offset + i * width for i in range(curN)]
#         #curCol = "g"
#         for i in range(len(hatch)):
#             rects1.extend(ax1.bar(curInd[i], np.mean(curList[i]), width, color=curCol, hatch=hatch[i]))#, yerr=np.std(curList[i]), ecolor='k'))
#         autolabel(ax1,rects1[int(2*el):],np.mean(curList,axis=1),np.std(curList,axis=1))
#
#         curList = []
#         if el==3: curList.append(getSingleData(fld, "INCOMP" , "PT_S", fluidData))
#         curList.append(getSingleData(fld, "HEOS"   , "DT_S", fluidData))
#         curList.append(getSingleData(fld, "REFPROP", "DT_S", fluidData))
#         lab.extend(["S"])
#         offset += step
#         curN   = len(curList)
#         curInd = [ offset + i * width for i in range(curN)]
#         #curCol = "g"
#         for i in range(len(hatch)):
#             rects1.extend(ax1.bar(curInd[i], np.mean(curList[i]), width, color=curCol, hatch=hatch[i]))#, yerr=np.std(curList[i]), ecolor='k'))
#         autolabel(ax1,rects1[int(3*el):],np.mean(curList,axis=1),np.std(curList,axis=1))
#
#         curList = []
#         if el==3: curList.append(getSingleData(fld, "INCOMP" , "PT_HDS", fluidData))
#         curList.append(getSingleData(fld, "HEOS"   , "DT_HPS", fluidData))
#         curList.append(getSingleData(fld, "REFPROP", "DT_HPS", fluidData))
#         lab.extend(["all"])
#         offset += step
#         curN   = len(curList)
#         curInd = [ offset + i * width for i in range(curN)]
#         #curCol = "g"
#         for i in range(len(hatch)):
#             rects1.extend(ax1.bar(curInd[i], np.mean(curList[i]), width, color=curCol, hatch=hatch[i]))#, yerr=np.std(curList[i]), ecolor='k'))
#         autolabel(ax1,rects1[int(4*el):],np.mean(curList,axis=1),np.std(curList,axis=1))
#
#         curList = []
#         if el==3: curList.append(getSingleData(fld, "INCOMP" , "HP_D", fluidData))
#         curList.append(getSingleData(fld, "HEOS"   , "HP_D", fluidData))
#         curList.append(getSingleData(fld, "REFPROP", "HP_D", fluidData))
#         lab.extend(["D"])
#         offset += step
#         curN   = len(curList)
#         curInd = [ offset + i * width for i in range(curN)]
#         curCol = ccycle.next()
#         for i in range(len(hatch)):
#             rects2.extend(ax2.bar(curInd[i], np.mean(curList[i]), width, color=curCol, hatch=hatch[i]))#, yerr=np.std(curList[i]), ecolor='k'))
#         autolabel(ax2,rects2[0:],np.mean(curList,axis=1),np.std(curList,axis=1))
#
#         curList = []
#         if el==3: curList.append(getSingleData(fld, "INCOMP" , "HP_T", fluidData))
#         curList.append(getSingleData(fld, "HEOS"   , "HP_T", fluidData))
#         curList.append(getSingleData(fld, "REFPROP", "HP_T", fluidData))
#         lab.extend(["T"])
#         offset += step
#         curN   = len(curList)
#         curInd = [ offset + i * width for i in range(curN)]
#         #curCol = "r"
#         for i in range(len(hatch)):
#             rects2.extend(ax2.bar(curInd[i], np.mean(curList[i]), width, color=curCol, hatch=hatch[i]))#, yerr=np.std(curList[i]), ecolor='k'))
#         autolabel(ax2,rects2[int(1*el):],np.mean(curList,axis=1),np.std(curList,axis=1))
#
#         curList = []
#         if el==3: curList.append(getSingleData(fld, "INCOMP" , "HP_S", fluidData))
#         curList.append(getSingleData(fld, "HEOS"   , "HP_S", fluidData))
#         curList.append(getSingleData(fld, "REFPROP", "HP_S", fluidData))
#         lab.extend(["S"])
#         offset += step
#         curN   = len(curList)
#         curInd = [ offset + i * width for i in range(curN)]
#         #curCol = "r"
#         for i in range(len(hatch)):
#             rects2.extend(ax2.bar(curInd[i], np.mean(curList[i]), width, color=curCol, hatch=hatch[i]))#, yerr=np.std(curList[i]), ecolor='k'))
#         autolabel(ax2,rects2[int(2*el):],np.mean(curList,axis=1),np.std(curList,axis=1))
#
#         curList = []
#         if el==3: curList.append(getSingleData(fld, "INCOMP" , "HP_DST", fluidData))
#         curList.append(getSingleData(fld, "HEOS"   , "HP_DST", fluidData))
#         curList.append(getSingleData(fld, "REFPROP", "HP_DST", fluidData))
#         lab.extend(["all"])
#         offset += step
#         curN   = len(curList)
#         curInd = [ offset + i * width for i in range(curN)]
#         #curCol = "r"
#         for i in range(len(hatch)):
#             rects2.extend(ax2.bar(curInd[i], np.mean(curList[i]), width, color=curCol, hatch=hatch[i]))#, yerr=np.std(curList[i]), ecolor='k'))
#         autolabel(ax2,rects2[int(3*el):],np.mean(curList,axis=1),np.std(curList,axis=1))
#
#
#         ids = np.arange(len(lab))
#         # add some text for labels, title and axes ticks
#         #if backend=="INCOMP": ax1.set_ylabel(r'Time per $f(p,T)$ call (\si{\us})')
#         if el==3: ax1.set_ylabel(r'Time per \texttt{Tmax} call (\SI{0.01}{\us}) and'+"\n"+r'per $f(p,T)$ and $f(\rho,T)$ call (\si{\us})')
#         else: ax1.set_ylabel(r'Time per \texttt{Tmax} call (\SI{0.01}{\us})'+"\n"+r'and per $f(\rho,T)$ call (\si{\us})')
#         ax2.set_ylabel(r'Time per $f(h,p)$ call (\si{\us})')
#
#         ax1.set_xticks(ids)
#         ax1.set_xticklabels([r"\texttt{"+i+r"}" for i in lab], rotation=0)
#         ax1.set_xlim([ids.min()-2.5*width,ids.max()+2.5*width])
#
#         ax1.spines['top'].set_visible(False)
#         ax2.spines['top'].set_visible(False)
#         ax1.xaxis.set_ticks_position('bottom')
#         ax2.xaxis.set_ticks_position('bottom')
#
#         handles = []
#         if el==3:
#             for h in (rects1[0], rects1[4], rects2[2]):
#                 handles.append(copy.copy(h))
#                 handles[-1].set_facecolor('white')
#                 handles.append(copy.copy(h))
#                 handles[-1].set_hatch('')
#             labels = (r'$p,T$-fit', r'\texttt{Tmax}', r'CoolProp', r'explicit, $f(p|\rho,T)$', r'REFPROP', r'implicit, $f(h,p)$')
#         else:
#             for h in (rects1[0], rects1[2], rects2[1]):
#                 handles.append(copy.copy(h))
#                 handles[-1].set_facecolor('white')
#                 handles.append(copy.copy(h))
#                 handles[-1].set_hatch('')
#             labels = (r'', r'\texttt{Tmax}', r'CoolProp', r'explicit, $f(\rho,T)$', r'REFPROP', r'implicit, $f(h,p)$')
#             handles[0] = mpatches.Patch(visible=False)
#
#         if bp is not None:
#             bp.drawLegend(ax=ax1,
#               loc='upper center',
#               bbox_to_anchor=(0.5, 1.4),
#               ncol=3,
#               handles=handles,
#               labels=labels)
#         else:
#             ax1.legend(handles,labels,
#               loc='upper center',
#               bbox_to_anchor=(0.5, 1.4),
#               ncol=3)
#
#         fg.savefig(path.join(getFigureFolder(),"TimeComp-"+fld+".pdf"))
#         if bp is not None: bp.savepgf(path.join(getFigureFolder(),"TimeComp-"+fld+".pgf"),fg,repList)
#         plt.close()
#
#     except Exception as e:
#         print(e)
#         pass

#############################################################
# The second figure compares the backend for the full calculation
#############################################################
backendExp = []
backendImp = []
fluidLabel = []
hatchLst    = ["///"  ,"\\\\\\"]
backendsLst = ["HEOS","REFPROP"]
for fluidstr in fluids[:-1]:
    _,fld,_ = splitFluid(fluidstr)
    #if fld=="CO2": fluidLabel.append("\ce{CO2}")
    #else:
    fluidLabel.append(fld)
    for backend in backendsLst:
        if backend not in fluidData[fld]:
            backendExp.append([0])
            backendImp.append([0])
            continue
        # Backend exists in fluid data
        try:
            inp = getInpList(backend)
            outExp = getTimeKey(inp[0], getOutList(inp[0])[-1])
            outImp = getTimeKey(inp[1], getOutList(inp[1])[-1])
            backendExp.append(getSingleData(fld, backend, outExp, fluidData))
            backendImp.append(getSingleData(fld, backend, outImp, fluidData))
        except Exception as e:
            backendExp.append([0])
            backendImp.append([0])
            print(e)
            pass

# Data is prepared, we can plot now.
if bp is not None:
    bp.figure = None
    fg1 = bp.getFigure()
    bp2 = BasePlotter()
    fg2 = bp2.getFigure()
    ccycle = bp.getColorCycle(length=3)
else:
    fg1 = plt.figure()
    fg2 = plt.figure()
    ccycle = cycle(["b","g","r"])

fg1.set_size_inches( (fg1.get_size_inches()[0]*1, fg1.get_size_inches()[1]*0.75) )
fg2.set_size_inches( (fg2.get_size_inches()[0]*1, fg2.get_size_inches()[1]*0.75) )
ccycle.next() # No incomp
#
#ax1 = fg.add_subplot(111)
#ax2 = ax1.twinx()
ax1 = fg1.add_subplot(111)
ax2 = fg2.add_subplot(111)

#entries = int(len(backendExp)/len(fluidLabel))
#one standard deviation above and below the mean of the data
rects1 = []
rects2 = []
col1 = ccycle.next()
col2 = ccycle.next()

numFluids   = len(fluidLabel)
numBackends = len(backendsLst)

step   = 1
width  = step/(numBackends+1)  # the width of the bars
offset = -0.5*numBackends*width

for f in range(numFluids):
    ids = np.empty((numBackends,))
    for b in range(numBackends):
        i = f*numBackends+b
        ids[b] = offset + f*step + b*width
        rects1.extend(ax1.bar(ids[b], np.mean(backendExp[i]), width, color=col1, hatch=hatchLst[b], rasterized=False))#, yerr=np.std(curList[i]), ecolor='k'))
        rects2.extend(ax2.bar(ids[b], np.mean(backendImp[i]), width, color=col2, hatch=hatchLst[b], rasterized=False))#, yerr=np.std(curList[i]), ecolor='k'))

y_max = np.max(np.concatenate((np.ravel(ax1.get_ylim()),np.ravel(ax2.get_ylim())/10.0)))

ax1.set_ylim([0,y_max])
ax2.set_ylim([0,y_max*10.0])

for ax in [ax1,ax2]:
    ax.set_xticks(range(numFluids))
    ax.set_xticklabels(fluidLabel, rotation=25)
    ax.set_xlim([0.0-0.5*step,numFluids-1+0.5*step])
    ax.spines['top'].set_visible(False)
    ax.xaxis.set_ticks_position('bottom')

    if   ax==ax1: rects = rects1
    elif ax==ax2: rects = rects2
    handles = (rects[0], rects[1])
    labels  = (r'CoolProp', r'REFPROP')
    anchor  = (0.5, 1.2)

    if bp is not None:
        bp.drawLegend(ax=ax,
          loc='upper center',
          bbox_to_anchor=anchor,
          ncol=numBackends,
          handles=handles,
          labels=labels)
    else:
        ax.legend(handles,labels,
          loc='upper center',
          bbox_to_anchor=anchor,
          ncol=numBackends)

ax1.set_ylabel(r'Time per $f(\rho,T)$ call (us)')
ax2.set_ylabel(r'Time per $f(h,p)$ call (us)')
fg1.savefig(path.join(getFigureFolder(),"TimeComp-backends-exp.pdf"))
fg2.savefig(path.join(getFigureFolder(),"TimeComp-backends-imp.pdf"))
if bp is not None:
    ax1.set_ylabel(r'Time per $f(\rho,T)$ call (\si{\us})')
    ax2.set_ylabel(r'Time per $f(h,p)$ call (\si{\us})')
    mpl.rcParams['text.usetex'] = True
    fg1.savefig(path.join(getFigureFolder(),"TimeComp-backends-exp-tex.pdf"))
    fg2.savefig(path.join(getFigureFolder(),"TimeComp-backends-imp-tex.pdf"))
    mpl.rcParams['text.usetex'] = False
    bp.savepgf(path.join(getFigureFolder(),"TimeComp-backends-exp.pgf"),fg1,repList)
    bp.savepgf(path.join(getFigureFolder(),"TimeComp-backends-imp.pgf"),fg2,repList)
plt.close('all')

#############################################################
# The third figure is a heat map of the execution times in
# log p h diagram
#############################################################
for fluidstr in fluids:
    try:
        _,fld,_ = splitFluid(fluidstr)
        for backend in fluidData[fld]:
            propsfluid = "::".join([backend,fluidstr])

            if backend!="INCOMP":
                TP = {}
                points = max(int(xypoints/2),250)
                T_range_TP = np.linspace(PropsSI('Ttriple',"T",0,"D",0,propsfluid)+1,PropsSI(  'Tcrit',"T",0,"D",0,propsfluid)-0.1,points)
                T_TP = np.append(T_range_TP,T_range_TP[::-1])
                Q_TP = np.zeros_like(T_TP)
                Q_TP[points:] = 1
                points *= 2

                out = ["D","H","P","S"]
                res = PropsSI(out,"T",T_TP,"Q",Q_TP,propsfluid)
                D_TP = res[:,0]
                H_TP = res[:,1]
                P_TP = res[:,2]
                S_TP = res[:,3]

                mask = np.isfinite(D_TP)
                if np.sum(mask)<points:
                    warn("There were not enough valid entries in your result vector. Reducing the number of points from {0:d} to {1:d}.".format(points,np.sum(mask)))
                    points = np.sum(mask)

                TP["T"] = T_TP[mask]
                TP["D"] = D_TP[mask]
                TP["H"] = H_TP[mask]
                TP["P"] = P_TP[mask]
                TP["S"] = S_TP[mask]
                #saveNpzData(TP)
            else:
                TP = None

            state = getStateObj(propsfluid)
            if backend=="HEOS" and state.has_melting_line():
                p_melt = np.logspace(np.log10(state.melting_line(CoolProp.constants.iP_min, CoolProp.constants.iT, 0)),np.log10(state.melting_line(CoolProp.constants.iP_max, CoolProp.constants.iT, 0)),xypoints)
                #p_melt = p_range
                ML = dict(T=[],D=[],H=[],S=[],P=p_melt)
                for p in p_melt:
                    try:
                        ML["T"].append(state.melting_line(CoolProp.constants.iT, CoolProp.constants.iP, p))
                    except Exception as ve:
                        ML["T"].append(np.inf)
                res = PropsSI(["D","H","P","S","T"],"T",ML["T"],"P",ML["P"],propsfluid)
                ML["D"] = res[:,0]
                ML["H"] = res[:,1]
                ML["P"] = res[:,2]
                ML["S"] = res[:,3]
                ML["T"] = res[:,4]
                mask = np.isfinite(ML["T"])
                ML["P"] = ML["P"][mask]
                ML["T"] = ML["T"][mask]
                ML["D"] = ML["D"][mask]
                ML["H"] = ML["H"][mask]
                ML["S"] = ML["S"][mask]
            else:
                ML = None
            #ML = {}

            IP = {}
            p_range,T_range = getPTRanges(propsfluid)
            critProps = getCriticalProps(propsfluid)
            try:
                IP["T"] = T_range
                IP["P"] = np.zeros_like(T_range) + critProps["P"]
                res = PropsSI(["D","H"],"T",IP["T"],"P",IP["P"],propsfluid)
                IP["D"] = res[:,0]
                IP["H"] = res[:,1]
            except Exception as ve:
                IP = None

            IT = {}
            try:
                IT["P"] = p_range
                IT["T"] = np.zeros_like(p_range) + critProps["T"]
                res = PropsSI(["D","H"],"T",IT["T"],"P",IT["P"],propsfluid)
                IT["D"] = res[:,0]
                IT["H"] = res[:,1]
            except Exception as ve:
                IT = None

            ID = {}
            try:
                ID["T"] = T_range
                ID["D"] = np.zeros_like(p_range) + critProps["D"]
                res = PropsSI(["P","H"],"T",ID["T"],"D",ID["D"],propsfluid)
                ID["P"] = res[:,0]
                ID["H"] = res[:,1]
            except Exception as ve:
                ID = None

            IH = {}
            try:
                IH["P"] = p_range
                IH["H"] = np.zeros_like(p_range) + critProps["H"]
                res = PropsSI(["D","T"],"P",IH["P"],"H",IH["H"],propsfluid)
                IH["D"] = res[:,0]
                IH["T"] = res[:,1]
            except Exception as ve:
                IH = None

            IS = {}
            try:
                IS["P"] = p_range
                IS["S"] = np.zeros_like(p_range) + critProps["S"]
                res = PropsSI(["D","H","T"],"P",IS["P"],"S",IS["S"],propsfluid)
                IS["D"] = res[:,0]
                IS["H"] = res[:,1]
                IS["T"] = res[:,2]
            except Exception as ve:
                IS = None

            for I in [IP,IT,ID,IH,IS]:
                if I is not None:
                    mask = np.isfinite(I["D"]) & np.isfinite(I["H"])
                    if np.sum(mask)<20: I = None
                    else:
                        for k in I:
                            I[k] = I[k][mask]






            for inp in getInpList(backend):
                if bp is not None:
                    bp.figure = None
                    fg = bp.getFigure()
                else:
                    fg = plt.figure()

                kind = getTimeKey(inp, getOutList(inp)[-1])
                t_data = getSingleData(fld, backend, kind, fluidData)
                x_data = getSingleData(fld, backend, "H", fluidData)
                y_data = getSingleData(fld, backend, "P", fluidData)

                gs  = gridspec.GridSpec(1, 2, wspace=None, hspace=None, width_ratios=[10,1])
                ax1 = fg.add_subplot(gs[0,0], axisbg='Tan')
                ax1.set_yscale('log')
                #ax2 = ax1.twinx()

                minHP = np.min(t_data)
                maxHP = np.max(t_data)
                minIT = np.percentile(t_data,10)
                maxIT = np.percentile(t_data,90)
                difIT = np.log10(maxIT/minIT)*0.25

                print(kind,": {0:7.2f} to {1:7.2f}".format(minHP,maxHP))
                if kind == "DT":
                    if bp is not None:
                        cx1=bp.getColourMap(reverse=True)
                    else:
                        cx1   = mpl.cm.get_cmap('cubehelix_r')
                    minHP = minIT
                    maxHP = np.power(10,np.log10(maxIT)+difIT)
                    #minHP = np.power(10,np.log10(np.percentile(t_data,10)*1e6))
                    #maxHP = np.power(10,np.log10(np.percentile(t_data,90)*1e6)*1.10)
                    #maxHP = np.power(10,1.10*np.log10(maxHP))
                    #minHP = np.percentile(t_data,10)*1e6
                    #maxHP = np.percentile(t_data,99)*1e6
                    #print(kind,": {0:7.2f} to {1:7.2f}".format(minHP,maxHP))
                    #minHP = 100
                    #maxHP = 20000
                else:
                    if bp is not None:
                        cx1 = bp.getColourMap()
                    else:
                        cx1   = mpl.cm.get_cmap('cubehelix')
                    minHP = np.power(10,np.log10(minIT)-difIT)
                    maxHP = maxIT
                    #minHP = np.power(10,np.log10(np.percentile(t_data,10)*1e6)*0.90)
                    #maxHP = np.power(10,np.log10(np.percentile(t_data,90)*1e6))
                    #minHP = np.percentile(t_data,01)*1e6
                    #maxHP = np.percentile(t_data,90)*1e6
                    #print(kind,": {0:7.2f} to {1:7.2f}".format(minHP,maxHP))
                    #minHP = 100
                    #maxHP = 20000




                #cx1_r = reverse_colourmap(cx1)
                cNorm  = mpl.colors.LogNorm(vmin=minHP, vmax=maxHP)
                #cNorm  = mpl.colors.LogNorm(vmin=ceil(minHP/1e1)*1e1, vmax=floor(maxHP/1e2)*1e2)
                #cNorm  = mpl.colors.Normalize(vmin=round(minHP,-2), vmax=round(maxHP,-2))
                colourSettings  = dict(c=t_data, edgecolors='none', cmap=cx1, norm=cNorm)
                pointSettings   = dict(s=6)
                scatterSettings = dict(rasterized=True, alpha=0.5)
                #scatterSettings = dict(rasterized=False, alpha=0.5)
                scatterSettings.update(colourSettings)
                scatterSettings.update(pointSettings)
                SC = ax1.scatter(x_data/1e6,y_data/1e5, **scatterSettings)

                for I in [TP,ML]:
                    if I is not None:
                        ax1.plot(I["H"]/1e6,I["P"]/1e5,lw=1.5,c='k')

                for I in [IP,IT,ID,IS,IH]:
                    if I is not None:
                        ax1.plot(I["H"]/1e6,I["P"]/1e5,lw=1.0,c='k',alpha=1)


                #ax1.set_xlim([0e+0,6e1])
                #ax1.set_ylim([5e-1,2e4])
                ax1.set_xlim([np.percentile(x_data/1e6,0.1),np.percentile(x_data/1e6,99.9)])
                ax1.set_ylim([np.percentile(y_data/1e5,0.1),np.percentile(y_data/1e5,99.9)])

                formatter = ticker.LogFormatter(base=10.0, labelOnlyBase=False)
                #formatter = ticker.ScalarFormatter()
                #ticks     = roundList(np.logspace(np.log10(ax1.get_ylim()[0]), np.log10(ax1.get_ylim()[1]), 5))
                #locator   = ticker.FixedLocator(ticks)
                #ax1.yaxis.set_major_locator(locator)
                ax1.yaxis.set_major_formatter(formatter)

                cax = fg.add_subplot(gs[0,1])
                formatter = ticker.ScalarFormatter()
                CB  = fg.colorbar(SC, cax=cax,format=formatter)
                CB.set_alpha(1)
                CB.locator = ticker.MaxNLocator(nbins=7)
                #ticks = roundList(np.logspace(np.log10(minHP), np.log10(maxHP), 5))
                #CB.locator = ticker.FixedLocator(ticks)
                CB.update_ticks()
                CB.draw_all()

                #fg.suptitle("f("+inp+")-"+backend.lower()+"-"+fld.lower())
                CB.set_label(backend.upper()+"::"+fld+', execution time per f('+inp[0]+","+inp[1]+') call (us)')
                ax1.set_xlabel(r'Specific enthalpy (MJ/kg)')
                ax1.set_ylabel(r'Pressure (bar)')
                fg.tight_layout()
                fg.savefig(path.join(getFigureFolder(),"TimeComp-"+inp+"-"+backend.lower()+"-"+fld.lower()+".pdf"))

                #CB.set_label(r'Execution time per call (\si{\us})')
                #ax1.set_xlabel(r'Specific enthalpy (\si{\mega\J\per\kg})')
                #ax1.set_ylabel(r'Pressure (\si{\bar})')
                #fg.tight_layout()
                #bp.savepgf(path.join(getFigureFolder(),"TimeComp-"+inp+"-"+backend.lower()+"-"+fld.lower()+".pgf"),fg,repList)
                plt.close()

    except Exception as e:
        print(e)
        pass
    plt.close('all')




















